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A granular media lattice gas (GMLG) model is used to 
study avalanches in a two-dimensional granular pile. We 
demonstrate the efficiency of the algorithm by showing that 
several features of the non-critical behaviour of real sandpile 
surfaces, such as the bounded outflow statistics or the finite- 
size effect of the time evolution of the pile mass, can be re- 
produced by this simulation approach. 



The concept of self-organized criticality (SOC) |l]] 
has become popular in the research of non-equilibrium 
systems. According to this idea, weakly driven non- 
equilibrium systems can spontaneously organize them- 
selves into a state of diverging characteristic length and 
time scales. It was suggested that avalanches on sandpile 
surfaces may fit into the SOC framework and this idea 
inspired considerable activity in the field of the physics 
of granular materials. 

Sandpiles soon proved to behave differently from the 
proposed picture in many aspects. Avalanches in a rotat- 
ing drum were found to occur at well-defined interwalls 
and the probability density of avalanche durations turned 
out to be sharply peaked Using a different experi- 
mental setup, where grains were dropped onto the top of 
a conical pile, Held et al. did find scaling in small piles 
||. However, if the size of the pile was above a certain 
value, the small fluctuations disappeared from the sys- 
tem. The crossover to this quasi-regular behaviour can 
be explained by the fact, that the surface of the pile has 
two characteristic angles rather than a single one: the 
angle of repose and the angle of marginal stability. In 
the case of sufficiently small piles the size of one single 
particle is large enough to raise the local angle above the 
critical value and therefore no hysteresis can be found 
Further measurements pointed out that the inter- 
vals between two consecutive large events contain several 
small avalanches with a power law size distribution ||. 
There have been attempts to predict the large events, 
but large avalanches appear to show the characteristics 
of a Markov-process [||. In accordance with this result, 
the power spectrum computed from the time sequence 
of avalanches follows a l// 2 rule (see e.g. By re- 

analyzing the results of several experimental investiga- 
tions, J. Feder has found that in the case of the small 
events stretched exponential functions (having a char- 
acteristic size) fit the avalanche size distributions more 



precisely. 

Recently two model systems have been found which do 
show SOC: in rice piles (experimental results: simula- 
tion: J9j) and avalanches of small particles in a bidisperse 
filling of a rotating drum [|l0| (simulation) . In both cases 
it turns out that criticality is a consequence of surface 
roughness, due to the elongated grains in the former case 
and the surface niches formed by the larger particles in 
the latter. It is also surprising that rice piles made up 
of grains of different aspect ratio - so varying only in a 
material constant in terms of critical behaviour - belong 
to different universality classes ||. 

The aim of this work is to present the Granular Me- 
dia Lattice Gas (GMLG) model through a study of 
avalanches in a pile. Experimental results often lack large 
enough statistics, molecular dynamics is inadequate for 
the same reason, while simplified computer models may 
exclude important aspects of real sandpile surfaces (such 
as the difference between the two critical angles or inertia 
effects). The simulation method we use here is a promis- 
ing compromise between the two: It has been designed to 
describe many features of real sandpiles while it is based 
on a fast parallel algorithm. 

The GMLG model is the generalization of a success- 
ful, fully discrete hydrodynamic algorithm pT|^2[ . It is 
defined on a triangular lattice, where indistinguishable 
point-like particles either travel at unit velocity along the 
lattice bonds or they can be at rest at the nodes and on 
the bonds. An update consists of the collision and prop- 
agation step. During the collision step the particles are 
scattered at the lattices nodes while during propagation 
step the moving particles are transferred to the nearest 
neighbour sites. 

While bulk collisions of the original hydrodynamic 
model conserve mass, energy and momentum, in the 
GMLG model energy dissipation and friction effects are 
also taken into account. As a result of the restricted set 
of velocities (0 or I), material parameters can be intro- 
duced only in a stochastic way, by means of probability 
variables. 

A driven sandpile is an example of such a granular 
system, where both fluidized and static regions can be 
present at the same time. The different behaviour of the 
material within these regions is reflected in the collision 
rules of the model as follows. 

In fluidized regions momentum is conserved, but col- 
lisions can dissipate energy. This is described by the 
parameter pk, which is the model equivalent of the en- 
ergy restitution coefficient. If energy is conserved at a 
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particular site - with a probability of pk - then the colli- 
sion rules are similar to that of the hydrodynamic model, 
with the exception, that rest particles may block possible 
scattering directions. If this is the case, we choose the 
after-collision configuration which provides the maximum 
mixing of states. If the collision is dissipative, then the 
maximum energy is dissipated while conserving momen- 
tum. We present three examples on Fig. |l| demonstrating 
how these principles work in the fluidized region. 

The compact static part of the pile behaves like a solid 
with a large mass, where friction effects are to be taken 
into account. When moving particles interact with the 
bulk, their momentum can be transferred through the 
force chains to the walls of the vessel. This is modelled 
by the friction variable p^. It gives the probability, with 
which a moving particle stops when arriving at a bulk site 
(see Fig. ||) . Bulk particles are such rest particles which 
are supported either by another bulk site or the bottom of 
the vessel. Although this definition itself does not guar- 
antee that there cannot be nodes that temporarily are 
not supported by the bottom or the walls - in principle, 
it would be possible to set up an algorithm that finds the 
bulk sites at the cost of efficiency - but in case of simple 
geometries our definition works correctly. (This process 
is similar to the capture process, introduced in a different 
context p|l.) 

The implementation of gravity for moving particles is 
straightforward, like for the hydrodynamic models [Q. 
Gravity rules applied to rest particles involve an effective 
metastability, in that moving particles can set off rest 
ones. This metastability is responsible for the hysteresis 
of the pile surface. Note, that in our case the "micro- 
scopic" rules result in the two different critical angles, 
we do not include them, what is the usual approach with 
cellular automaton sandpile models (e.g. [|ip,|i~4|] ) . 

After applying the collision, friction and gravity rules, 
the particles are transferred to the nearest neighbour 
sites. The GMLG propagation step differs from the hy- 
drodynamic model in one aspect. Dissipative binary col- 
lisions may take place here, since there can be up to two 
particles on each bond between NN sites, as it is illus- 
trated on Fig. ||. 

With the model described above and with similar ap- 
proaches it is possible to simulate such different systems 
as pipe flows, shaken boxes, granular mixtures or static 
piles 15-18|]. In this paper we examine the properties of 
avalanches. 

As a first test of the model we calculate the depen- 
dence of the angle of repose of the pile on the friction 
variable p^. We have chosen a method described in [fl9|| , 
the steady-state filling of a silo and the result is shown 
on Fig. The error bars give information about the 
difference between the angle of repose and the angle of 
marginal stability. The angles are rather high compared 
to real- world systems, which is a consequence of the un- 
derlying lattice (the smallest angle of repose is 30° in the 
model) . 

The system studied - a 2D box with one open side - 



is a classical setup for examining avalanches and is also 
analogous with the Hele-Shaw cell used e.g. by Frette et 
al ||. The particles are dropped near to the side- wall in 
such a way that a half-pile is building up (see Fig. ||). 
The pile is driven quasi-statically, that is the particles are 
added one by one, after all activity caused by dropping 
the previous particle ceased. The size L of the system 
is defined by the length of the horizontal support. The 
measurements start after the stationary state has been 
reached. 

First we consider the time evolution of the total mass 
of the pile (see Fig. ^|). The time unit here is the interval 
between dropping two consecutive grains, which can last 
from one single update to several hundred updates long. 

The graphs on Fig. ^| represent two runs, where all 
material parameters are kept identical, while the system 
sizes are varied. Although the sizes are modified only by a 
factor of two, the curves are qualitatively different. The 
first one (L = 24) is a function irregularly fluctuating 
on many time scales, but the second one (L = 48) is 
much more regular, quasi-periodic with a period of about 
4000 timesteps. This can be underlined by comparing the 
power spectra of the two data series (see Fig. |?]). In case 
of the larger pile a peak develops, which corresponds to 
a frequency of 1/4000. This finite-size effect is in nice 
agreement with the experimental findings of Held et al 
S. This result also makes it possible to calibrate the 
length of the model system to experimental scales. 

As a next step we study the distributions of avalanches. 
We use two quantities which are sufficient to characterize 
the size of an avalanche: the T lifetime of an event (in 
update units) and the number of particles falling off the 
support, that is the M mass of a droplet. Note that M 
does not contain information about the small avalanches 
not reaching the rim of the support. A probability den- 
sity curve contains data obtained from typically 10 6 -10 7 
updates. In (20) it was found for a SOC automaton that 
the outflow statistics has multiscaling properties. The 
quality of our generated data was not good enough to 
check for this property. Therefore, in a simple finite-size 
scaling framework we look for the distributions of these 
quantities in the following form: 



p(T,L)~ L-ofriT/L* 



p(M,L)~L-Pf M (M/L») 

Fig. U and Fig. [)] show the probability densities of 
these quantities for different system sizes. All curves 
were logarithmically binned and rescaled using the ansatz 
above, but for comparison the inset on Fig. || displays 
two typical raw data curves. 

The first apparent feature is that the finite-size ef- 
fect observed in the time evolution of the pile appears 
in the avalanche statistics too. Both the lifetime and 
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the droplet probability densities have a power law form 
(straight line on a log-log plot) with a sharp cutoff for 
small system sizes (L < 40) while a pronounced peak de- 
velops for larger sizes. This behaviour was also observed 
in the IBM experiment ||] . This feature is somewhat less 
obvious in the statistics of the duration times, since the 
sharper characteristic peaks are smoothed when binning 
the curves, but they are well-defined, see the original data 
curves on the inset of Fig. ||. 

The good data collapse demonstrates, that the finite- 
size scaling ansatz is well suited for the duration time 
distributions. A well-defined exponent can be found for 
the small avalanches: p(T, L) ~ T y with y = 1.92 ± 0.05. 
The numerically found exponents are also consistent with 
the y = —a/A criterion, which follows from the properties 
of the /r(a) scaling function for a — > 0. Note, that the 
value of y is very close to the critical exponent found in 
a so-called local version of a class of rice pile models || . 
The droplet distributions, however, are different in the 
two models. 

The outflow statistics of the small, L = 20 pile needs 
more attention. Although the distribution of small 
droplets can be described by a power law, it is also possi- 
ble to fit the whole curve, including large avalanches, with 
a stretched exponential ansatz of the form f(M, 20) ~ 
exp(— c * M 1 ). The exponent (7 = 0.34) that was used 
by Feder |?J while re-analyzing experimental data gives 
an excellent fit (see Fig. |l0|) . 

Finally we compare the GMLG model to a recent ex- 
periment concerning rice piles The friction rule in 
the simulation can be interpreted as the probability for a 
grain to get trapped in a local surface minimum, there- 
fore a higher p^ (in terms of the rice pile experiment a 
higher p^ can be regarded as a higher grain aspect ratio) 
should result in a rougher surface and this is what can ac- 
tually be seen in our model. This may suggest a smooth 
transition from the sandpile-type distribution to that of 
the rice pile. By increasing the coefficient of friction, 
the outflow probability distribution does get broader, but 
there is no change in the scaling properties for different 
Pfi values. Higher p^ means a larger angle of repose and 
also a larger difference between the critical angles (Fig. 
||), therefore the big avalanches consist of more particles. 
Fig. [n] shows outflow distributions for two different fric- 
tion values. 

Although the droplet and lifetime distributions are suf- 
ficient to describe the avalanche statistics, for the sake 
of an easier comparison we have also checked the dis- 
tribution of the dissipated energy in an avalanche. By 
comparing the total potential energy of the pile before 
and after an avalanche the dissipated energy can be ob- 
tained. The distributions for four system sizes are plot- 
ted on Fig. 12. The finite-size effects found in lifetime- 
and droplet distributions can be observed here, too. The 
power law part of the probability densities have an expo- 
nent of 1.63 ±0.06, which is different from the one found 
experimentally for rice piles (« 2.0). The model expo- 
nent is very close to the value measured in an another 



cellular automaton rice-pile model 0. 

In summary, the GMLG-pile has many features in com- 
mon with real sandpiles. The non-trivial finite-size ef- 
fects, such as the characteristically different time depen- 
dence of the pile mass at different sizes, the outflow statis- 
tics are in good qualitative agreement with experimen- 
tal findings. With a proper calibration of the geometry 
the model may even give quantitatively good results. It 
is important to emphasize that the GMLG-model does 
not involve any built-in assumptions about the nature 
of avalanches (as opposed to many cellular automaton 
models). 
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FIG.0 Three examples of the collision rules in fluidized 
regions. The left column shows three pre-collision con- 
figurations. The middle and the right column show con- 
figurations - together with their probabilities - after an 
clastic and an inelastic collision, respectively. Cases a), 
b) and c) are examples of a head-on collision between 
two particles, with zero, one and two rest particles also 
present on the site. In case a), if the collision is elastic, 
one of two possible configurations are chosen with equal 
(l/2pfe) probability. With 1—pk probability the collision 
is inelastic and the particles stop. In case b) one scatter- 
ing direction is blocked, which is free in example a). In 
case c) both directions are blocked and binary collisions 
between particles on opposite bonds take place. 

FIG-D An example for the application of the friction 
rule, when a moving particle arrives at a bulk site. (The 
marked particles belong to the bulk.) With probability 

the particle looses its energy and becomes part of the 
bulk. In the opposite case its momentum is conserved 
and the node is not considered to belong to the static 
part any more. 

FIG.|| Illustration of the propagation step. Since two 
particles can be present on a bond, binary collisions may 
take place. Those particles after an inelastic collision 
are marked. The figure shows examples for all possible 
configurations on the bonds. 

FIG. [| Dependence of the angle of repose of the pile on 
the friction parameter. 

FIG. ^ Simulation snapshot of the pile. Note, that each 
circle represents a lattice node, which can be occupied by 
up to seven particles. 

FIG-D The total mass (number of particles) of two piles 
vs. time. (One timestep is an interval long enough to 
contain the avalanches of the longest lifetime, rather than 
an update unit). (The friction parameter is p^ — 0.34). 
A factor of two in the system size results in significantly 
different time sequence. 

FIG. [7] The power spectra computed from total mass 
- time graphs for two different sizes. (Both curves have 
been averaged for 25 runs.) The peak at / = 2.5 10 -4 ( 
At » 4000) in case of L — 48 shows the appearance of a 
characteristic time scale. 

FIG. H Finite size scaling of lifetime distributions (p p = 
0.5). The best collapse is obtained at the scaling expo- 
nents A = 0.8 ± 0.02, a = 1.6 ± 0.04. The exponent of 
the power law part of the distribution is y = 1.92 ± 0.05. 
All distribution curves on this graph and on the following 
ones are binned. 

FIG. I Finite size scaling of droplet distributions (p^ = 
0.5). The j3 = 1.4, /i = 0.7 exponents here used for visu- 
alisation only, since it is apparent that a scaling according 
to the ansatz in the text cannot be applied. 



FIG .|| A streched exponential fit f(M, 20) ~ exp(-c* 
Af 7 ) for the droplet probability density, in the case of 
a small pile. 7 = 0.34 is used here, which is equal to 
the fit used for experimental data. (The value of c is 
2.0 ± 0.2 in the model, which depends on the details of 
the dynamics.) 

FIG.pd|The effect the friction coefficient on the droplet 
distribution. The system size is L = 80. 

FIG.p"2| Probability density of the dissipated energy in 
an avalanche for four system sizes. The exponent of the 
power-law part of the curves is 1.63 ± 0.06. 
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